%这个octave脚本简易演示了波穿过界面时的透射与反射，使用了基于显式Euler法的PDE数值求解
%CC-BY 4

clc
clear

n=1000;
dx = 0.01;
dt = 0.001;
c=[1;2]; %wave speed
T=1; %period
w=2*pi/T;

u0 = zeros(n,1);
u1 = zeros(n,1);
u2 = zeros(n,1);

k=0;
imgindex = 0;

Areflection = (c(2)-c(1))/(c(1)+c(2)); %基于波动力学理论计算的反射系数
Atransmission = 2*c(2)/(c(1)+c(2)); %透射系数

_x = linspace(0,dx*(n-1),n);
_x = _x';

figure()
for t = 0:dt:10
##  for i = 2:n-1
##    if i < n/2
##      j=1;
##    else
##      j=2;
##    endif
##    u2(i) = 2*u1(i) - u0(i) + (c(j)*dt/dx).^2*(u1(i-1)-2*u1(i) + u1(i+1));
##  endfor

  % 向量化
  mp = ceil(n/2); %middle_point
  u2(2:mp-1) = 2*u1(2:mp-1) - u0(2:mp-1) + (c(1)*dt/dx).^2*(u1(1:mp-2)-2*u1(2:mp-1) + u1(3:mp));
  u2(mp:n-1) = 2*u1(mp:n-1) - u0(mp:n-1) + (c(2)*dt/dx).^2*(u1(mp-1:n-2)-2*u1(mp:n-1) + u1(mp+1:n));

  if t > T
    u2(1) = 0;
  else
    u2(1) = sin(w*t);
  endif
  u2(n) = 0;

  u0 = u1;
  u1 = u2;

  if mod(k,20) == 1
    clf
    hold on
    axis equal
    axis([0 n*dx -2 2])

    plot(_x,u2)
    line([0 n/2*dx],  [Areflection Areflection],'linestyle','--')
    line([n/2*dx n*dx],  [Atransmission Atransmission],'linestyle','--')
    line([n/2*dx n/2*dx], [-2 2]) %interface

    pause(0.01)
    k;
##    print(['wave_interface_' num2str(imgindex) '.png'],'-dpng');
    imgindex++;
  end
  k++;
end
